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I.  DESIGN  PHILOSPHY 


IPROP  is  a  nonlinear,  three-dimensional,  relativistic,  electromagnetic 
particle-in-cell  (PIC)  beam  simulation  code  written  to  treat  high  current 
electron  beam  propagation  in  the  atmosphere.  As  such.  It  has  all  the  usual 
features  of  a  plasma  simulation  code  plus  air  chemistry  and  particle  scat¬ 
tering.  The  IPROP  field  solver  has  been  modified  to  accommodate  air  con¬ 
ductivity,  large  aspect  ratio  spatial  zoning,  and  a  moving  coordinate 
mesh.  Some  specialized  diagnostics  have  been  added  as  well.  The  field 
solver,  chemistry  routine,  particle  dynamics  package,  and  diagnostics  pro¬ 
cedures  are  described  in  Section  II  through  V,  respectively.  Here,  we  com¬ 
ment  briefly  on  the  rationale  for  the  particular  design  chosen  for  IPROP. 


At  the  time  that  IPROP  was  written.  In  early  1983,  the  few  existing 
nonaxl symmetric  propagation  codes  were  limited  to  simulating  the  m=l  hose 
Instability  in  a  single  transverse  plane^*^CL  Higher  azimuthal  modes, 
arising  from,  for  example,  fl lamentation,  were  not  treated.  Moreover,  the 
m=l  particle  dynamics  and  air  chemistry  were  linearized.  The  usual  frozen 
field  and  paraxial  particle  motion  approximations  also  were  employed^ 
For  a  variety  of  reasons  these  approximations  and  limitations  were  Incom¬ 
patible  with  our  needs. 


The  development  of  IPROP  had  been  commissioned  by  Sandla  National  Lab¬ 
oratories  as  part  of  the  RADLAC  high  current  beam  propagation  program.  The 
RADLAC  beam  was  expected  to  be  annular  and,  therefore,  subject  to  rapid 
fl lamentation.  VISHNU5  (and,  more  recently,  IBEX6)  beam  experiments 
indicated  that  fl lamentation  typically  does  not  disrupt  propagation  but  can 
cause  the  beam  radial  profile  to  fill  in  partially.  As  a  consequence, 
IPROP  had  to  be  capable  of  following  arbitrary  azimuthal  modes  nonl In- 
early.  That  the  RADLAC  beam  would  rotate  required  that  dynamics  In  both 
transverse  planes  be  kept.  The  frozen  field  and  paraxial  particle  motion 
approximations  were  dropped  to  permit  simulations  of  near  term  experiments 
involving  beam  energies  of  a  few  to  several  MeV. 


It  should  also  be  mentioned  that  the  PIC  linearized  hose  codes 
suffered  from  noise  and  numerical  Instability  problems  at  that  time.  One 
significant  source  of  numerical  noise  Is  easily  Identified.  Even  In  the 
absence  of  applied  forces,  perturbed  particle  positions  Increase  linearly. 
The  effect  of  this  secular  growth  on  the  perturbed  fields  vanishes  In  the 
limit  of  an  Infinite  number  of  particles,  but  not  uniformly.  As  time 
increases  In  a  computation,  progressi vely  more  particles  are  required  to 
achieve  a  given  level  of  cancellation  among  the  secular  contributions  of 
the  perturbed  particles  positions  to  the  perturbed  currents.  Equivalently, 
for  a  fixed  number  of  particles  the  noise  level  grows  with  time.  This 
apparently  Is  the  cause  of  the  late-time  failure  of  phase-mix  damping  some¬ 
times  observed  In  such  codes.  (This  difficulty  can  be  overcome  by  rein¬ 
itializing  perturbed  particle  quantities  every  so  often,  as  Is  done  In 
SIW1.3)  Nonlinear  simulation  codes  certainly  are  not  Immune  to  numer¬ 
ical  noise.  However,  years  of  experience  with  PIC  codes  In  a  variety  of 
applications  Indicate  that  noise  usually  Is  not  a  severe  problem  when  rea¬ 
sonable  precautions  are  taken.  It  may  be  that  nonlinear  codes  are  quieter 
at  late  times,  because  the  greatest  spatial  separation  that  two  particles 
Initially  nearby  In  phase-space  can  attain  Is  the  diameter  of  a  particle 
orbit. 

The  occasional  exponential  growth  of  perturbed  particle  orbits  In 
self-pinched  equilibria  even  without  perturbed  fields,  reported  in 
RINGBEARER  II  simulations,7  Is  less  easily  explained.  A  simple  calcula¬ 
tion  shows  that  such  behavior  Is  nonphysical.  We  have  speculated  that  the 
runaway  perturbed  positions  may  be  due  to  stochastic  numerical  resonance 
hopping,8  but  we  have  not  pursued  this  Idea.  In  any  event,  it  did  not 
seem  prudent  to  develop  IPROP  as  a  linearized  PIC  code  In  the  face  of  this 
problem  of  unknown  origin. 

Practical  considerations  also  Influenced  the  design  of  IPROP.  Ear¬ 
lier,  we  had  written  the  two-dimensional  beam  propagation  code  CPR0P.9 
which  successfully  employed  algorithms  not  Involving  the  frozen  field  or 


paraxial  particle  approximations.  CPROP  was  developed  quickly  and  at  mini¬ 
mal  cost  from  the  general  purpose  two-dimensional  plasma  simulation  code 
CCUBE.10  IVORY,11  the  three-dimensional  extension  of  CCUBE,  was, 
therefore,  the  natural  basis  for  IPROP.  Basically,  creating  IPROP  required 
transferring  the  CCUBE  modifications  to  IVORY  and  then  devising  a  method  of 
handling  the  nonlinear  coupling  among  azimuthal  modes  of  the  conductivity 
and  electromagnetic  fields.  The  latter  was,  of  course,  not  trivial. 
Upgraded  diagnostics  were  added  only  recently. 

II.  THREE-DIMENSIONAL  FIELD  SOLVER 


The  CPROP  axlsymmetrlc  field  solver  was  designed  to  be  numerically 
accurate  and  stable  for  arbitrary  scaler  conductivities,  have  an  unre- 
strlctlve  Courant  time  step  limit,  allow  computations  In  both  the  labora¬ 
tory  and  beam  frames,  and  accommodate  low  as  well  as  high  energy  particle 
beams.  Its  nonaxl symmetric  generalization  Is  meant,  as  well,  to  treat  any 
modestly  sized  set  of  not  necessarily  consecutive  azimuthal  modes  In  either 
or  both  transverse  planes.  Following  the  axlsymmetrlc  derivation,12  we 
cast  Maxwell's  equations  Into  coupled  forward-  and  backward-going  wave 
equations. 


i[it  +  (1  “  v)  !d  <Ee  -  Br>  -  utf  -  (1  +  v)  m  (Ee  +  V 
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The  axial  mesh  velocity  v  is  assumed  non-negative  for  coding  convenience. 

The  six  equations  are  next  partially  Integrated  along  their  axial 
characteristics  and  finite  differenced  axially  and  radially.  We  obtain 

OT. 

f1-')]  *  C'll  -  %)  -  f  M)  lrh  B"z+1>1 

■  £r  (E  -  “2)  (E?  -  BS)  *  Hr  H2  (E?,m  -  B"e,m) 


5 


■r1-1  (i 

2  2 x2f  + 

*%(> 

-  «2)  bz 

«i)  b;-1 

where 

•Cj  *  min 

(¥-«) 

(l-v)At 

Az 


1  B_  Fn+l,1 
r  Br  r  fce 


At 


I  cn+1»1 
r  B0  fcr 


At 


and  the  Index  palp  p  Is  given  by  (n+1,1+1)  or  (n,i)  depending  on  whether 
W2  assumes  its  first  or  second  value  in  Eq.  (13).  In  the  equations  At  is 
the  time  step  and  Az  Is  the  axial  cell  size.  The  Indices  n  and  1  Indicate 
time  step  and  axial  cell  numbers,  respectively,  with  1  Increasing  from  the 
tail  of  the  beam  to  the  head.  See  Fig.  1. 


Centered  radial  differences  are  represented  by  differentials  In  Eqs. 
( 7 ) — ( 12)  for  notational  simplicity.  Er,  Be,  and  Bz  are  located  radially  at 
cell  edges,  while  Br,  Eq,  and  Ez  are  at  cell  centers.  The  first  three 
fields  are  required  to  vanish  on  axis.  The  second  three  are  set  to  zero  at 
the  outer  radial  edge  of  the  mesh.  For  simulation  of  open  air  propagation, 
the  outer  radius  Is  chosen  to  be  at  least  twenty  times  the  nominal  beam 
Bennett  radius.  Radial  zoning  uniform  In  the  coordinate  £, 

$  *  a  An(l  +  r/a)  (15) 


with  a  the  Bennett  radius,  normally  is  employed  to  yield  cells  uniform  In 
size  within  the  beam  and  expanding  linearly  at  larger  radii. 

The  azimuthal  dependences  of  the  field  and  current  components  and  of 
functions  of  the  conductivity  are  given  by 
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Figure  1.  Typical  z-t  coordinate  meshes  for  (a)  Az  <  (l+v)At,  and  (b)  Az  > 
(l+v)At.  Forward-  and  backward-going  light  lines,  interpolation 
quantities  Wi  and  W2,  and  forward  light  line  integration  time 
also  are  shown.  Omitted  are  the  conductivity  characteristic  and 
integration  time  ti. 
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in  Eqs.  (7)-(12).  Consistent  with  this  finite  Fourier  expansion,  products 
of  f(o)  and  the  fields  or  currents  are  replaced  by  convolutions,  (Treat¬ 
ment  of  the  convolutions  is  discussed  in  the  Appendix.)  We  emphasize  that 
only  azimuthal  modes  of  interest  need  be  kept  In  a  particular  computation. 
For  Instance,  simulations  of  filamentation  In  the  recent  IBEX  low  density 
air  propagation  experiments®, 13  might  employ  only  m  *  0,  4,  8. 

Equations  (7)-(12)  with  auxiliary  conditions  (13)-(18)  are  solved 
numerically  at  each  time  step  by  sweeping  axially  from  the  head  of  the  beam 
back  (1.  e.,  decreasing  1).  At  each  axial  slice  the  right  sides  of  the 
equations  are  evaluated  from  quantities  already  known  from  either  the  pre¬ 
ceding  time  step  or  the  axial  disk  next  further  forward  in  the  beam,  as 
Indicated  in  Fig.  1.  Er  and  Br  can  be  eliminated  from  the  left  sides 
of  the  equations,  leaving  4»Nr»Ne  coupled  equations  and  a  like  number 
of  unknowns;  l.e.,  of  order  103  equations  and  unknowns  for  typical  num¬ 
bers  of  radial  zones  Nr  and  azimuthal  modes  Ne.  Although  inverting 
such  a  system  Is  possible,  doing  so  would  be  both  slow  and  memory  inten¬ 
sive.  Instead,  azimuthal  modes  are  decoupled  by  treating  the  convolutions 
iteratively,  allowing  the  resulting  systems  of  4»Nr  equations  for  each 
mode  to  be  Inverted  as  in  CPROP.  Beginning  with  m  =  0,  we  evaluate  higher 
m  field  contributions  to  the  convolutions  using  values  from  the  preceding 
tlmestep.  Ez,  Be,  Bz,  and  Ee  for  m  =  0  then  are  computed  from  the  radial 
finite  difference  equations  by  Gauss  elimination,  and  the  corresponding 
Er  and  Br  are  determined  by  back  substitution.  For  the  next  higher 
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azimuthal  mode,  we  employ  the  new  field  values  for  m  *  0  but  the  old  values 
for  the  remaining  modes.  The  radial  equations  then  are  Inverted  just  as 
for  m  *  0.  This  procedure  Is  repeated  for  the  remaining  azimuthal  modes.  In 
each  case  evaluating  the  convolutions  with  the  most  recent  field  values 
available.  One  may  sweep  through  the  modes  as  often  as  necessary  to  obtain 
the  desired  accuracy.  In  practice,  we  find  that  a  single  pass  is  suffi¬ 
cient.  This  procedure  is  straightforward  and  reasonably  fast. 

Had  the  convolutions  not  been  treated  iteratively,  the  Courant  condi¬ 
tion  associated  with  Eqs.  ( 7 )- ( 12)  would  be 

At  <  Az/(l-v)  (19) 

Probably,  the  Iteration  procedure  makes  the  actual  Courant  condition  more 
restrictive,  but  this  has  caused  no  apparent  difficulties.  For  v  near 
unity.  At  Is  constrained  by  particle  dynamics,  as  in  the  frozen  field 
approximation. 

Fields  are  initialized  at  the  beginning  of  a  simulation  based  on  the 
frozen  field  approximation. 

III.  AIR  CHEMISTRY  ROUTINE 


I PROP,  when  first  written,  employed  the  PHOENIX  air  chemistry  algo 
rithm  developed  by  Lawrence  Livermore  National  Laboratory.^ 


do 

dt 


K  0  +  v  a 


2 

a  a 


(19) 


In  our  dimensionless  system  of  units,  the  impact  Ionization,  avalanche  ion¬ 
ization,  and  recombination  coefficients  are  given  by 

K  =  0.593  (20) 


v  = 


_ 1.163*105  PS3 _ 

1  +  2.667* 101  S  +  2.242*105  S*  +  6.916*102  S3 


(21) 
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a  * 


(22) 


3.2«10~4  P 
(1  +  5.23* 105  S)0,7 

with  S  =  (E/P)2,  the  electric  field  E  scaled  to  511  kV/cm,  and  the  air 
density  P  measured  in  atmospheres  (at  standard  temperature  and  pressure). 

In  the  course  of  delta-ray  studies  described  in  Sec.  2.4,  we  found  it 
useful  to  replace  the  constant  impact  ionization  coefficient  K  in  Eq.  (20) 
by  a  more  accurate,  energy  dependent  expression  approximating  Bethe's  elec¬ 
tron  energy  loss  formula. 15 

K  =  0.302  +  0.054  Jin  y  (23) 

where  y  is  the  beam  electron  average  energy,  redetermined  periodically. 
(Using  the  local  rather  than  the  average  energy  would  be  slightly  more 
accurate  but  somewhat  more  expensive  computationally.)  Equation  (23) 
reduces  to  Eq.  (20)  at  about  200  MeV.  At  lower  energies  the  new  expression 
has  a  slight  stabilizing  effect  on  the  beam  by  enhancing  impact  ionization 
at  the  expense  of  avalanche.  For  the  low  current,  low  energy  VISHNU  exper¬ 
iments,  Eq.  (23)  was  necessary  even  to  obtain  correct  return  current  decay 
times.8 


Recently,  the  recombination  coefficient  a  was  replaced  with 


_ 1.19«10"4  P 

(1  +  35.4  S1/4  +  44.0  s1/2)0*39 


Although  we  have  not  yet  had  the  opportunity  to  assess  the  effects  of  this 
change,  it  is  clear  that  Eq.  (24)  leads  to  faster  recombination  for  larger 
values  of  E/P  and,  therefore,  should  be  stabilizing.  At  the  same  time,  a 
temperature  dependent  momentum  transfer  cross  section  was  introduced  by  the 
simple  procedure  of  multiplying  the  conductivity  by 
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r _ h02 _ 

(1  +  35.4  S1/4  +  44.0  s1/Z)0,36 

efore  using  it  in  the  field  solver.  Eq.  (24)  and  (25)  have  been  cali- 
, rated  against  the  HICHEM  code  at  high  air  densities  for  a  wide  range  of 
ieam  currents. 


The  air  conductivity  is  computed  at  half-integer  timesteps  so  that  it 
:an  be  advanced  with  the  fields  in  a  leap-frog  fashion.  The  computation  is 
lone  in  real  space  to  avoid  convolutions  and  then  decomposed  into  azimuthal 
nodes  for  use  in  the  field  solver.  The  finite  difference  equation  corres- 
jonding  to  Eq.  (19)  is 


a0 

i  +  R+v  1 

L1  +  w  \ 

Pi 

1  +  5r  1 

(•N 

+  ( 

.*> ') 

a0  =  (l  -  W3j  a*  +  W3  on,1+1 

R2  =  v2  +  4  a  K  |j| 
with  tj  defined  in  Eq.  (13)  and 


(26) 


(27) 

(28) 
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min 


vAt  \ 
AZ  ) 


(29) 


.  .  (27),  the  index  pair  \x  is  given  by  (n+l,i+l)  or  (n,i)  depending  on 

whett..  W3  assumes  its  first  or  second  value.  Note  that  Eq.  (26)  is 
exact  for  constant  fields. 


IV.  PARTICLE  DYNAMICS 


The  basic  particle  transport  algorithms  are  essentially  identical  to 
those  used  in  IVORY*!  and  so  will  not  be  discussed  here.  Bilinear  inter¬ 
polation  is  employed  in  the  particle-field  interface. 
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RADIAL  DIFFERENCES  HANDLED  IN  STRAIGHTFORWARD 
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PARTICULAR  FORM  OF  CONDUCTIVITY  INTEGRATING 
FACTOR  CHOSEN  TO  REPRODUCE  CORRECT  RESISTIVE 
DECAY  RATES 
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FIELD  ALGORITHM  DERIVED  FROM  MAXWELLS  S 
EQUATIONS  CAST  IN  TERMS  OF  FORWARD, 
BACKWARD  GOING  WAVES 
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AZIMUTHALLY  VARYING  FIELDS  TREATED  FLEXIBLY  BY 
pamtitw ATlftN  OF  FOURIER.  REAL  SPACE  ALGORITHMS 
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(QT  •  F  •  Q)  (QT  •  E) 


(A-5) 


Hence,  multiplying  the  transformed  fields  by  the  conductivity  tensor 
qT«F»q  generates  the  desired  convolution.  Note  that,  although  this  dis¬ 
cussion  is  framed  in  terms  of  a  spatial  mesh,  the  mesh  Is  not  actually 
required  for  most  of  the  manipulations. 


Examples  of  similarity  transformation  matrices  are 
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;  M  =  2 
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APPENDIX 


The  azimuthal  finite  Fourier  expansions  and  resulting  convolutions 
associated  with  Eq.  (16)-(18)  are  treated  explicitly  as  follows.  Suppose 
that  the  electromagnetic  fields  E  and  B,  and  the  currents  J  are  defined  on 
a  uniformly  spaced  azimuthal  mesh  of  M  cells.  The  fields  satisfy  differ¬ 
ential  equations  containing  the  operator 


D 


(A-l) 


The  finite  difference  representation  of  (A-l)  can  be  diagonalized  by  means 
of  a  similarity  transformation  on  Its  eigenvectors,^  evaluated  on  the 
mesh  and  normalized. 


(A-2) 


The  unnormalized  eigenvectors  are 


,  ,  (cos  (n-1)  mA0  m  >  0 

*<">  =1  1  m  =  0 

r"  (  -sin  (n-1)  mA0  m  <  0 


(A-3) 


Here,  A0  =  0max/M*  Applying  the  transformation  to  D«E  yields 

(QT  •  D  •  Q)  (QT  •  E)  =  <-mZ>  (QT  •  E)  (A-4) 

where  qT»e  Is  the  transformed  field  and  <-m2>  is  a  diagonal  matrix  with 
elements  -m2. 


Products  of  the  form  f(c)«E,  appearing  in  Eq.  (7),  (9),  and  (11)  are 
handled  In  an  analogous  fashion.  Let  F  be  a  diagonal  matrix  consisting  of 
f(o)  evaluated  at  the  M  azimuthal  mesh  points.  Employing  the  similarity 
transformation  described  above  to  F»E  gives 


including  comparisons  with  experiment,  has  been  described  In  greater  detail 
22 

elsewhere.  Scattering  Is  Implemented  In  IPROP  by  applying  a  deflec¬ 
tion  to  the  beam  electrons  every  few  cm  of  propagation  in  full  density  air 
and  proportionately  less  often  at  lower  densities.  The  deflection  angle  is 
chosen  randomly  from  a  large  set  of  previously  computed  small  angles  form¬ 
ing  a  truncated  Mollere  distribution  or  from  an  analytical  expression  for 
occasional  large  angles. 

V.  DIAGNOSTICS 

Satisfactorily  displaying  three-dimensional  phenomena  in  two-dimen¬ 
sional  plots  Is  Intrinsically  difficult.  IPROP  diagnostics  are  for  the 
most  part  simple  generalizations  of  CPROP  diagnostics.  Improved  output  is 
a  subject  of  continuing  research. 

Detailed  particle  behavior  Is  depicted  by  various  two-dimensional  pro¬ 
jections  In  six-dimensional  phase  space.  Energy  and  momentum  histograms 
and  time  histories  of  average  momenta  also  are  provided.  Individual  azi¬ 
muthal  modes  of  the  electromagnetic  fields  and  of  the  beam  current,  net 
current,  and  conductivity  distributions  are  displayed  In  contour  and  cross 
section  plots.  Time  histories  of  all  these  quantities  at  specified  loca¬ 
tions  on  the  computational  mesh  are  available  together  with  their  Fourier 
spectra.  The  amplitude  and  position  of  the  electric  field  spike  also  are 
traced  in  time.  Various  radial  moments  of  the  beam  and  net  current  have 
been  added  recently.  Other  recent  additions  to  the  diagnostics  are  des¬ 
cribed  In  Ref.  23. 

Output  usually  takes  the  form  of  microfiche  and  16  mm  color  movies. 
Paper  plots  and  35  mm  color  slides  also  are  available. 


Several  hundred  beam  electrons  are  loaded  Into  each  axial  slice  of  the 
coordinate  mesh  during  Initialization  to  give  a  Bennett  radial  pro- 
fl le16,17 


j2  *  -7  (» + <3°> 

na  '  a  r 

with  an  axially  dependent  half-radius  a.  Electron  axial  momenta  are  picked 
from  the  corresponding  relatlvistlcally  Invariant  distribution 

P2/t 

f(pz)  ■  T  -pTT -  p„  >  Pz  >  0  <3I> 

e  0  -  1 


Specifying  the  beam  energy  then  determines  the  transverse  momenta.  It  Is 
easy  to  show  for  high  energy  beams  that  the  resulting  transverse  momentum 
distribution  Is  approximately  Maxwellian  with  the  RMS  average  (p0T)l/2. 

The  mechanics  of  randomly  selecting  momenta  for  Eq.  (31)  Is  simplified  by 
the  usual  trick  of  equating  f(pz)  dpz  to  dC,  where  C  Is  a  set  of  random 
numbers  uniformly  distributed  between  zero  and  one.  There  results 


*  T  X n  [l  +  C  (eP°/T-l)] 


(32) 


Electrons  can  be  loaded  pairwise  to  minimize  statistical  fluctuations  at 
early  times.  The  electrons  of  each  pair  are  separated  azimuthally  by  n  but 
otherwise  have  Identical  phase  space  coordinates. 

Electron  scattering  by  air  molecules  is  modeled  using  the  Moll  ere 
formal  Ism, 18  which  consists  basically  of  Williams  scattering!?  at  small 
angles  and  Coulomb  scattering  at  large. 20  As  compared  with  the  usual 
Rossl-Grelsen  formalism, 21  Mollere  scattering  leads  to  a  smaller  radial 
expansion  rate  combined  with  a  very  slow  loss  of  beam  particles.  The 
Impact  of  various  scattering  models  on  computed  Nordsieck  lengths. 
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FIELD  EXPANDED  AZIMUTHALLY  IN  SINES,  COSINES 


M  a  •  Ma  <■ 


FIELD  EQUATIONS  INVERTED  BY  COMBINATION  OF 
EXPLICIT,  ITERATIVE,  AND  IMPLICIT  METHODS 
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)P  EMPLOYES  E/P  AIR  CONDUCTIVITY  MODEL 
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AIR  CHEMISTRY  NUMERICAL  ALGORITHM  EXACTLY 
SATISFIES  MODEL  EQUATION  FOR  CONSTANT  SOURCES 


CONDUCTIVITY,  FIELDS  EVALUATED  ALTERNATELY  IN  TIME 


PARTICLE  EQUATIONS  SOLVED  WITHOUT  LIMITING 
APPROXIMATIONS 
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ORBIT  AVERAGING  DAMPS  NUMERICAL  FLUCTUATIONS 
REDUCES  COMPUTER  COSTS 
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COMPREHENSIVE  DIAGNOSTICS  STILL  A  SUBJECT 
OF  CONTINUING  RESEARCH 
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IPROP  VERSATILITY  DEMONSTRATED  IN  DIVERSE 
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BEAM  TRANSPORT  IN  IFR  CHANNELS  ALSO  STUDIED 
WITH  IPROP 
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SIGNIFICANT  IPROP  ENHANCEMENTS  PLANNED 
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